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ABSTRACT 

This study presents a method for approximating the multidimensional effects of Rayleigh-Taylor 
instability as a modification of the one-dimensional hydro equations. This modification is similar to 
the Shakura-Sunyaev a prescription for modeling the coarse-grained effects of turbulence in astro- 
physical disks. The model introduces several dimensionless tunable parameters that are calibrated 
by comparing with high-resolution two-dimensional axisymmetric numerical calculations of Rayleigh- 
Taylor unstable flows. A complete description of the model is presented, along with a handful of 
test problems that demonstrate the extent to which the one-dimensional model is able to reproduce 
multidimensional effects. 

Subject headings: hydrodynamics — turbulence — shock waves — instabilities — ISM: jets and out¬ 
flows — ISM: supernova remnants — supernovae: general 


1. INTRODUCTION 

Many a strophysical outflows are Rayleigh-Taylor (RT) 
unstable (jChevalier fc Kleiiil 119781 1. This instability af¬ 
fects the properties of the flow, significantly enough that 
one-dimensional (ID) calculations assuming spherical 
symmetry may not match up with direct observations of 
supernova ejecta. RT causes disruption of sharp density 
jumps at contact discontinuities, mixing of ejecta with 
the circumstellar medium (CSM), interactions between 
the unstable region and the reverse shock, and line broad¬ 
ening due to turbulent fluctuations. This turbulence 
may also genera te magnet i c fields via small-scale tur- 
bulent dynamo l)Jun et al.l Il995t iDuffell fc MacFadvenl 
[Soil), and might even alter the forward shock propa- 
gation, if there is signihcant cooling in the shock fr ont 
(jBlondin fc Ellisonll2001l : IDuffell fc MacFadvenll20lJl . 

Because of the importance of these multidimen¬ 
sional effects, many numerical investigations have been 
launched over the past several decades to determine 
the properties of RT-unstable flows. The very first 
two-dimensional (2D) n umerical calculations were by 
iChevalier fc Kleiiil (jl978l l. who studied structure forma¬ 
tion in the context of Type II supernovae. Higher- 
resolution 2D results became possible decades later, 
and the problem w as more accurately tackled by 
IChevalier et al.l (|I992ll , who ca lculated gr o wth ra tes both 
analytically and numerically. Linn et ah! (jI995ll studied 
two and three-dimensional (3D) RT with magnetic fields 
in an idealized context, finding that magnetic fields af¬ 
fected the g rowth rate and that RT easily amplihes mag¬ 
netic fields. lJun h Normanl (|I996ll demonstrated in the 
supernova context (in 2D) that RT could cause these 
fields to align with turbulent struct ures, affecting po - 
larization of synchrotron emission. iKane et aH (j200flf l 
studied the difference between 2D and 3D, but still in a 
local sense (looking at single-mode perturbations). That 
study found that the growth of RT is 30 — 35% stronger 
i n 3D than in 2D. 

iBlondin fc Ellison! ()2001h first demonstrated the im¬ 
portance of cooling on the dynamics of RT. Specifi- 
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cally, cosmic rays provide signihcant cooling in shocks, 
resulting in shallower pressure gradients, allowing the 
Rayleigh-Taylor hngers to catch up to the forward shock. 
This was demonstrated hrst by varying the adiabatic 
index, showing a dramatic change in the dynamics. 
The importance of cosmic rays o n the dynam i cs ha s 
also been show n obse rvatioi iallv (iWarren et al.l 120051) . 
iFraschetti et al.l l) 201011 and iFerrand et al.l ( ~010ll uer- 
formed the hrst 3D global studies of RT, that also in¬ 
cluded a prescribed model for cosmic ray cooling (rather 
than varying the adiabatic index). RT has also been 


studied in the relativistic case, in the context of g£ 
rav bursts ([Duffell & MacFadvenll2013l 1201411. 

On the experimental side, RT has been v 
explored by various groups (e.g. IBudil et al.l 

imma 

widely 

199a 

Remington et al.lll997lll999l:IDimonte & Schneider! 

2000t 

Robev et al.l 12001b ISmalvuk et al.l 12005b IKuranz et al. 

2010b iCasner et al.ll2012ll. This is partly due to its im- 


portance in inertial conhn e ment f usion (|Lindl fc MeadI 
119751: iVerdon et al.l 119821: iLindll 1199511 . Richtmyer- 
Meshkov instability is also important in thi s context 
([Clendinning et ani2003l : iThornber et al.ll2011ll . 

Civen RT’s importance, and given that many numer¬ 
ical studies of supernova dynamics are still carried out 
in ID, one would expect the development of a ID model 
that could be used by theorists numerically studying su¬ 
pernovae using ID radiation hydrodynamics codes. Sur¬ 
prisingly, the only method which could be easily intro¬ 
duced into a standard I D hydro code was developed over 
40 years ago (|Cnlllll97^ . before multidimensional numer¬ 
ical studies of RT were possible. 


There have been numerous studies which attempt to 
theoretical ly model the nonlinear effects o f RT in various 


wavs 

'e.ff. iHaanl 119891: lAlon et alJ 11994 iDimonte et al. 

2004 

; 

), 

Shvarts et al.l 11995b lOfer et al.l Il99a lOron et al. 

2001 

using both experimental and numerical results. 


More general turbulence models have been suggested, 
e. g. the “One-D imensional Turbulence” formulation 
of iKersteinl (Il999f) (see also iWunsch fc Kerstehil 120011 : 
lAshurst fc Kersteinll2005ll , which uses a Monte Carlo ap¬ 
proach in an attempt to model detailed statistical prop¬ 
erties of the turbulence, including a Kolmogorov 5/3 cas- 
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cade. Most of these models generate a series of ordinary 
differential equations (ODEs) describing the evolution of 
various properties of a mixing layer. In contrast, Gull’s 
model was expressible as a set of partial differential equa¬ 
tions (PDEs) (adding additional terms to Euler’s equa¬ 
tions), which is useful for incorporation into existing ID 
hydro codes. However, only a handful of PDE turbulence 

models have been built specifically for KD _ 

A few notable such RT models include lYounesI (1198911 . 
a multiphase model for experiments using m ultiple im- 
miscible fluids, and the ODE formulation of Ramshaw 

S , which has also been cast in PDE form ( Ramshaw 
Both of these models have two drawbacks which 
make it difficult to apply them to the supernova con¬ 
text. First, both model the system using a multi¬ 
phase approach, suggesting an initial condition of two 
separate fluids, and hence foreknowledge of the po¬ 
sition of the unstable layer. Ideally, an RT model 
should allow unstable layers to grow anywhere they 
happen to emerge, treating every fluid element in the 
system equally. Such an approach requires a single¬ 
fluid model. Secondly, an expanding outflow sets lim¬ 
its on the maximal allowed coherent length scale in 
the turbulence, and this does not ap pear to be ac ¬ 
co unted for i n eith er of the models of lYoungsl l)1989f) 
or iRamsha^ l)2000ll . Kerstein’s ODT model was also 
applied to both RT and Richtmeyer-Meshkov instability 
(jGonzalez-.Iuez et al.ll201^ Llozefik et al.l[20Ififl . but this 
model has similar drawbacks as the models of Youngs 
and Ramshaw. Furthermore, the version applied to RT 
requires specific knowledge of the position of the mix¬ 
ing layer, and tracks its growth with time as part of the 
model; it is not clear how one would extend this to a 
model where multiple mixing layers could spontaneously 
form anywhere that instability is triggered. Moreover, 
the full reproduction of all details of the statistics of the 
turbulence may be overkill for the purposes of a ID su¬ 
pernova remnant calculation. 

In practice, when evolving supernovae in ID, one is of¬ 
ten forced to artificially “smear out” regions which are 
RT unstable in post-processing fe.g. iPint o fc Wo oslevI 
Il988t iKasen fc Woosle\l I2009I iHeger fc WooslevI l201(m . 
Such an operation is obviously undesirable, but it may 
be necessary if nothing else is done to model the effects 
of RT. 

Building on Gull’s 1973 method, an improved ID 
model is presented in this work, one that is informed and 
calibrated by multidimensional numerical RT calcula¬ 
tions. Following Gull’s original idea, an additional scalar 
quantity is evolved along with the usual hydro variables. 
This field, k, represents the magnitude of turbulent 
RT fluctuations, analogous to the Shakura-Sunyaev a 
prescription in disk physics (jShakura fc Sunvaevlll973ll . 
This is also a similar approach to that employed to deal 
with convective regions in ID stellar evolution codes, 
and turbulence-induced m ixing in type la supernovae 
(IWooslev et al.l[2n09l[20lTll . 

K is prescribed a growth rate in RT-unstable regions of 
the flow, and is used to calculate a local diffusion con¬ 
stant, 77 , that causes mixing of all conserved variables. 
This surprisingly simple prescription is enough to pro¬ 
duce dynamics reasonably consistent with true multidi¬ 
mensional RT calculations, sufficiently so that the model 
could be considered an improvement to a purely ID cal- 



Fig. 1. — Cartoon depiction of the unstable contact discontinuity 
between ejecta and CSM in the outflow from a supernova. 

culation. 

This one-dimensional model is detailed in Section [2l 
and is calibrated and tested in several supernova contexts 
in Section [3l A summary is presented in Section |4l 

2. MODEL DESCRIPTION 

This study will produce an augmented system of hydro 
equations, in an attempt to reproduce some of the coarse¬ 
grained effects in a one-dimensional model. The unmod¬ 
ified spherical ID hydro equations (in conservation-law 


form) are 

dtip)+ r~‘^ir‘^ipv)y = 0 ( 1 ) 

dt{pv) + + P))'= 2P/r (2) 

dtietot) + + P)v)' = 0, (3) 

where primes denote radial derivatives, p is density, v 
is the radial velocity, P is pressure, and ttot is the total 
energy density, defined as: 

etot = + f-Th (4) 

where tTh is the thermal energy density, and the equa¬ 
tions are closed by an equation of state: 

p = (^-l)eTh, (5) 


where 7 is the adiabatic index (7 = 5/3 here). Addi¬ 
tionally, one can evolve passive scalars (e.g. representing 
various nuclear abundances) as an additional conserva¬ 
tion law: 

dt{pX)+r~‘^{r‘^{pXv))'= Q. (6) 

In this study, a passive scalar X will be used to mark 
the separation between ejecta and the CSM. X=1 for 
the ejecta and X=0 for the CSM. This will be used to 
measure mixing between the two fluids. 

The RT instability occurs at large density gradients, 
typically at contact discontinuities between two fluids. 
In the supernova context, the two fluids are typically 
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ejecta (mass ejected in the explosion) and the GSM (ex¬ 
ternal mass swept up by the outflow). As the ejecta 
sweeps up the CSM, shocks are generated by the colli¬ 
sion. A forward shock rushes ahead into the CSM, and a 
reverse shock propagates back into the ejecta, notifying 
fluid elements in the ejecta t hat it is time to decelerate 
(|McKeelll974lChevalieilll982h . In the hot region between 
the two shocks resides the contact discontinuity, separat¬ 
ing ejecta from CSM (Figure [J). For typical outflows 
at early times, this contact discontinuity is unstable (see 
also Figure [2]). 

The instability is driven by the effective gravitational 
force felt in the decelerating reference frame of the con¬ 
tact discontinuity. In practical terms, unstable configu¬ 
rations exist whenever the pressur e gradient has an oppo - 
site sign to the density gradient (jChandrasekharlflQBlh . 
In this case, the pressure gradient is positive, as the 
ejecta is decelerating, and at early times the ejecta is 
more dense than the CSM, so the density gradient is 
sharply negative across the contact discontinuity, at least 
until a significant amount of the CSM is swept up. 

In order to account for these multidimensional effects, 
an additional scalar k is evolved with the ID flow, k 
represents the relative strength of turbulent fluctuations, 
similar to the Shakura-Sunyaev a viscosity. 

The scalar quantity n will be interpreted as 


and therefore turbu lent decay will no t be accurately cap¬ 
tured in this study. iMac Lo-(^ (|1999D showed that turbu¬ 
lence damps efficiently, and therefore this is an important 
consideration. 

The goal is to augment the ID equations m with ad- 
ditonal terms which mimic the effects of turbulence. Al¬ 
though the turbulence is subsonic, it propagates through 
steep density gradients, so that there are not just velocity 
perturbations Su, but also density fluctuations Sp. These 
fluctuations are significant enough that they are the dom¬ 
inant source of diffusion of all quantities (as opposed to 
the standard Reynolds stress terms p{SuiSuj) <C P). 

For example, consider the equation of mass continuity, 

dtp + V ■{pv) = 0. (10) 

the mass flux pv can be rewritten in terms of mean and 
fluctuating quantities: 

Fp = {p'^ = (Pavgi4,vg + <5pWavg + PavgSu + 5 p5u) (11) 

where the subscript “avg” denotes averaged quantities 
(this subscript will now be dropped for simplicity). In 
the overall average, the first term is the bulk mass flux, 
and the second and third terms average to zero. The last 
term 


« = (7) 

where 5u represents turbulent fluctuations of velocity 
about the background mean (averaged over a sufficiently 
large patch), and Cs is the local sound speed, = jP/p- 
Note that k is of order the ratio of kinetic to thermal 
energy in the turbulence: 


6Fp = (SpSu) (12) 

is nonzero if density fluctuations are correlated with ve¬ 
locity fluctuations. They should indeed correlate, if there 
is a density gradient in the flow. Assume that the fluc¬ 
tuation in density is given by the density that the fluid 
element would have carried one correlation time ago: 


K = pSu^/{jP) 


2 ^pSu^ 
7 P 


2 sk 
7(7 - 1) CTh 


( 8 ) 


In stable regions of the flow, k will be evolved as a 
passive scalar: 


dtipn) + r ^{r^{pKv)y = 0. (9) 

Note that this implicitly assumes that k stays fixed, 
even as a fluid element expands or contracts. This as¬ 
sumption is based on the idea that random kinetic fluc¬ 
tuations should have an effective adiabatic index of 5/3, 
and should therefore maintain a constant proportion with 
the thermal energy. 

This also neglects many source terms which would ap¬ 
pear in a formal Reynolds decomposition of the fluid 
equations. Viscous dissipation is neglected, as this en¬ 
tire treatment assumes zero viscosity, due to the large 
Reynolds numbers of these systems. Shear production 
and Buoyant production of turbulence are also neglected 
here, and these terms should be included in a more com¬ 
plete model, which would more accurately describe the 
evolution of the turbulence. In the present study, it is 
assumed that Rayleigh-Taylor instability dominates the 
growth of turbulence. 

Equation (O also does not take into account decay due 
to turbulent dissipation. This can potentially be added 
as a source term, as will be described later. It bears men¬ 
tioning, however, that the present work is restricted to 
2D numerical calculations for calibration of this model. 


Sp^—rSu-Vp, (13) 

where r is a correlation time. Given this assumption, 
mass flux can be re-expressed as 

Fp = pv — T (Su^) p' = pv — rjp'. (14) 

where p = t {Su^') is a diffusion constant. Similar ar¬ 
guments can be made to build diffusive fluxes for all con¬ 
served quantitiefl To summarize, neglecting growth or 
decay of turbulence, the equations should now have the 
form 


dt{p)+r '^{r'^{pv - pp'))' = 0 (15) 

dt{pv) + r~^{r‘^{pv‘^ + P — p^pv)'))' = 2P/r (16) 

5t(etot) + ?'"2(r2((etot + P)v - pe'tot))' = 0. (17) 

dtipX) + r-^r^ipXv - p{pXy)y = 0 (18) 

dt{pK) + r~'^{r^{pKV — p^pk)'))' = 0 (19) 


^ It should be noted that the precise form of these terms de¬ 
pends on assumptions made on how the conserved densities change 
as they are carried along by a fluid element. In this study, such 
effects were neglected, but future improvements should include 
these effects; for example, as a subsonic fluid element propagates 
through a pressure gradient, the density adjusts to ma intain spe¬ 
cific entropy. Therefore, a more accurate form for m should be 
-rPV75^. V(pP-i/7) 
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where rj = t (^6iP) = A(|(5m|) = X^/kCs, where A is a 
correlation length, discussed below. 

Growth of K is assumed to only take place in unsta¬ 
ble regions. Instability occurs wherever the gradients of 
pressure and density have opposite signs. It should be 
noted that k is only meant to represent the fluctuations 
of the largest-scale modes, as these modes dominate the 
mixing and turbulent stresses. The length scale A of these 
largest modes can be estimated by the following exercise: 

In order to have coherent structures in an expanding 
flow, the flow must not be expanding too quickly for the 
turbulent velocities to merge structures together on a 
given scale. The largest-scale modes will be at the critical 
scale such that the relative expansion velocity between 
two fluid elements separated by A in the background flow 
is of order the turbulent velocities: 

AV • u ~ (5 m (20) 

Assuming the background flow is expanding radially at 
a velocity of order the sound speed, and velocity fluctu¬ 
ations are given by 

V-?;~Cs/r, (5 m ~-v/k Cs. (21) 

This provides an estimate for the largest coherent scale 
of the turbulence: 

A ~ i/k r. (22) 

Note that this estimate may be inaccurate for a number 
of reasons (for example, the incorrect assumption that 
the flow expands at the sound speed, when actually there 
can be a large disparit y between sound speeds across t he 
contact discontinuity (|Truelove &: McK^ 119991 l2000fl l. 
In reality, any estimate of A from local quantities is prob¬ 
ably flawed, as A can cover large scales, and is not nec- 
essarily local l y calculable. Previous ID RT models (e.g . 
lYoungsiri989l : lRamsha-^l200r)l : iGonzalez-.Inez 6111111201311 

were designed in an idealized “infinite slab” context, 
where arbitrarily large wavelengths are allowed. Here, 
the spherical geometry and the expanding flow set the 
scale A. 

The nonlinear growth rate of Rayleigh-Taylor can be 
determined by imagining a mixing layer with density 
jump between two fluids with density pi and p 2 - Under 
constant acceleration g, the width of the mixing layer, 
denoted by h, grows quadratically with time in the non¬ 
linear regime (jDimonte et al.l [2004( 1: 

h{t) = agAt^ (23) 

where A= {pi — p 2 )l(,Pi + P 2 ) is the Atwood number, 
g is the acceleration felt by the mixing layer, and (a is a 
constant, empirically determined to have a value between 
0.03 and 0.08. The growth rate of this mixing layer is 
then 

Trt = h/h = 2/t = 2^JagA/h. (24) 

Assume that the density varies linearly from pi to p 2 
over the width h. Given that A/h = Ap/{2ph), (where 
Ap is the density jump and p = (pi -I- P2)/2 is the aver¬ 
age density in the mixing layer), and Ap/h = —p', the 
derivative of density in the mixing layer, and g = P'/p, 
this reduces to a growth rate of 


Trt - V-P'P'/P^- (25) 

where the constants have been dropped, and the minus 
sign is necessary as RT only grows when p' and P' have 
opposite signs. This growth rate translates into a source 
term for Equation ([5]): 


S+ = {A + B^i)rRTP = (A + Br) ^-P'p\ (26) 

where A and B are arbitrary dimensionless constants, 
and the source term is only added in regions where the 
pressure and density gradients have opposite signs, i.e. in 
regions where the square root evaluates to a real number. 
Note that this gives exponential growth for sufficiently 
large k, > A/H, but the small-amplitude growth phase is 
linear. This is because the goal is to model the nonlin¬ 
ear growth phase of the large-scale turbulence in super¬ 
nova remnants, after a significant amount of GSM mass 
has been swept up by the ejecta. In numerical studies 
it has been found that the satura ted instability is inde¬ 
pend ent of the seed perturbation (|Dnffell k-. MacFadvenI 
nnn), and therefore small-scale turbulence amounts to 
some nonzero growth rate at early times, after which ex¬ 
ponential growth takes over. This might not apply at 
early times in the evolution of the supernova. 

Turbulent decay can be taken to be exponential, with 
a decay timescale proportional to an eddy turn-over: 

Tdecay ^ A/|(5m|. (27) 

where A is the characteristic eddy size. Therefore, a 
sink term S~ should also be added: 

S~ = -Dpn/ Tdecay = -DpKCs/r. (28) 

in the current study, this term is turned off, D = 0, 
as the results are only compared with 2D turbulence, 
which should not exhibit significant decay. The term is 
included in the equations, however, so that future 3D 
studies can provide a decay constant. Finally, the value 
of K is used to calculate a kinematic diffusion constant 
rj. This diffusion constant appears in equation (1141) as 
rj T or 

Tj ^ \Su\X. (29) 

(this is also clear from dimensional grounds). Given A 
and using |(5m| = ^/kCs, 

g = C{KCsr), (30) 

where C is another dimensionless constant that will be 
calibrated using multidimensional calculations. 

All conserved quantities (including pn) are mixed ac¬ 
cording to the diffusion constant g. To summarize, the 
complete augmented ID system of equations can be ex¬ 
pressed in conservation-law form as 

dtU + r-^{r^{F-gU')y = S, (31) 

giving one equation for each conserved variable U: 

U = {p, pv, etot, pX, pk} . (32) 

The fluxes are given by 
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Fig. 2. — RT turbulence in a calculation of the idealized supernova at time t = 1 (Test #1, Section l3. II color denotes log of density). 


F = [pv, pv'^ -1- P, {etot + P)v, pXv, pnv} 

(33) 

and the source terms are 

,S={0,2P/r,0,0,5+ + 5-}, 

(34) 

where 

S'^ = (A-I-i?K)-\/max(—P'p', 0), 

(35) 

S~ = -Dpncs/r, 

(36) 

r] = CkCsT. 

(37) 


The choice of A is nontrivial, as it is an inherently 
global property of the flow, and therefor e cannot be com¬ 
puted from local quantities. iGu H (1197,11) effectively made 
the choice that A be associated with the zone size A ~ Ar, 
but this is of course problematic, as the solution might 
then depend on resolution. 

In this study, the assumption is made that A = ^/Kr, 
an assumption that is reasonably well-motivated based 
on the arguments above, but an assumption nonetheless. 
On the other hand, reasonable results are still possible 
for appropriate choices of the constants A, B, C, and D. 
A more accurate model may be possible using a better- 
informed prescription for A. For the choice A = y/l^r, 
two-dimensional RT is most accurately approximated 
by the constants A = 1.18 x 10“^, B = 1.2, C = 0.102, 
D = 0, as will be demonstrated in the next section. 

Implementation in an existing ID hydro code is very 
straightforward, especially if the code explicitly evolves 
the conservative form of the equations. For a given con¬ 
served quantity, U, the flux of U is simply modified by 
the replacement 

FiU) F{U) - TjU'. (38) 

All other aspects of the code (e.g. Riemann solvers) 
remain unchanged. This way of adding the diffusive term 


is explicit, which would be problematic for large diffusion 
constants, but as n is typically small ~ 1 — 10 %, this does 
not pose a major problem. The only additional issue is 
a correction to the timestep; in addition to the Courant 
condition, one must enforce the criterion 

At < Ar^/ry. (39) 

Again, for small 77 this is not a major problem. Such 
timestep restrictions might also be avoided using an im¬ 
plicit time integrator. 

3. NUMERICAL TESTS 

The numerical tests in this study are all performed in 
a similar way. Three calculations are performed for each 
test: a ID calculation with no RT model, a 2D axisym- 
metric calculation with real turbulence, and a ID calcu¬ 
lation with the RT model specified above. These tests are 
used for both calibration and evaluation of the method. 
Upon performing the first test, the following values for 
the tunable constants are found: A = 1.18 x 10“^, 
B = 1.2, C = 0.102, D = 0. These chosen values for the 
constants are used for all three numerical tests. It should 
be noted that the dynamics were found to be largely in¬ 
sensitive to the value of A, so long as it is small but 
nonzero. This could be linked to the fact that the large- 
scale dynamics of RT have been no ted to be largely insen¬ 
sitive to initial seed perturbations (|Duffell fc MacFadvenI 
nnn). However, this statement may depe nd on the wave¬ 
lengt h of the initial seed perturbation (|Dimonte et al.l 

The 2D axisymmetric calculati ons are performed us¬ 
ing the high ly accurate JET code (|Duffell fc MacFadvenI 
[ 2 OIII Holl). JET is a moving-mesh hydrodynamics 
code which is specifically tailored to the study of as- 
trophysical jets and outflows. The numerical method 
is based on the relativisti c moving-mesh TESS code 
l|Duffell fc MacFadvenll20lil) . but with a more restricted, 
radially shearing mesh. Computational zones can move 
and shear radially, so that each radial track behaves 
much like a one-dimensional Lagrangian code, and the 
different tracks are coupled laterally by transverse fluxes. 
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So me additional details on the technique are presented 
in iDuffell fc MacFadvenI l)2013f) . Mesh motion can sub¬ 
stantially reduce errors from advecting the turbulence 
across the computational domain. Another advantage of 
the moving-mesh technique is that the inner and outer 
boundaries follow the blast, so that the flow can move 
over a large dynamic range, while the code only needs 
to cover a modest dynamic range at any single time. 
Multidimensional calculations have a typical resolution 
of A9 ~ Ar/r Ri 4 x 10“'^. In the appendix, a resolu¬ 
tion study is performed demonstrating that this resolu¬ 
tion is more than sufficient to capture the coarse-grained 
properties of the turbulence. Mesh refinement and de¬ 
refinement is employed so that computational zones are 
typically kept at an aspect ratio close to unity. ID calcu¬ 
lations are performed using a custom-built ID moving- 
mesh hydro code. The ID source code is publicly avail¬ 
able at https://github.com/duffell/RTlD. 

The outer boundary condition is Dirichlet, setting the 
values at this boundary equal to the initial condition. 
This is appropriate as the outer boundary always con¬ 
sists of unshocked material, which is stationary until 
the ejecta collides with it. Mesh refinement and de¬ 
refinement is performed based on the aspect ratio of 
zones. If a zone has an aspect ratio higher than 5:1 
or lower than 1:5, the zone is either split in half, or 
merged with its neighbor. In either case, this action is 
performed in such a way as to conserve mass, energy, 
and momentum during the refinement or de-refinement 
process. In principle, it is possible to seed the instabil¬ 
ity with an initial perturbation. In practice, it has been 
found that the nonlinear phase is insensitive to the initial 
seed field, so long as this pe rturbation is small (see also 
IDuffell fc MacFadvenI l201,Tl) . This study works in the 
limiting case where the initial perturbation is arbitrarily 
small, and the instability is seeded by grid noise. The 
fact that the coarse-grained end state is converged (see 
Appendix) demonstrates that the final state is indeed 
insensitive to the magnitude of the initial perturbation. 

The first test, the “Idealized Supernova” is used to 
calibrate the ID model, determining appropriate choices 
of constants A, B, and C by comparing with a very well- 
resolved 2D calculation. This test will also be used to 
measure the effectiveness of the model, but these results 
should be taken with a grain of salt since the model was 
calibrated to perform as well as possible on this problem. 

The “Soft Explosion” is a test which is designed to pro¬ 
duce smoother, less dramatic shocks, to facilitate com¬ 
parison with numerical schemes which are not as robust 
at shock capturing. This initial condition is smooth ev¬ 
erywhere and the pressure is non-negligible so that the 
shocks generated are not as strong. The contact wave in 
this case is not a sharp discontinuity, and can be easily 
resolved. 

The “Sustained Engine” starts with a static initial con¬ 
dition, and deposits thermal energy on a non-negligible 
timescale. This test provides several independent checks 
on the model, including the fact that the engine gener¬ 
ates an RT-unstable region from the opposite direction 
(negative pressure gradient and positive density gradi¬ 
ent). This test also includes a parameter which controls 
the steepness of the contact discontinuity, which will be 
varied to further test the robustness of the model. 

The “Brief Encounter” is a somewhat more compli- 
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Fig. 3.— ID angle-averaged profiles of the idealized supernova. 
In this figure, density is plotted alongside the passive scalar X. This 
passive scalar is designed to measure mixing of ejecta with CSM 
(X = 1 for ejecta, and X = 0 in the CSM). 


cated initial condition, used to test the effectiveness of 
the model in less idealized scenarios. The ejecta den¬ 
sity in this model obeys a broken power-law in radius, 
instead of a purely uniform ejecta model, and the CSM 
is described by a shell of material that is only encoun¬ 
tered by the ejecta for a brief amount of time. This test 
and the final test are also calculated in cgs units, al¬ 
though it should be stressed that the source terms and 
diffusive fluxes were introduced in a scale-invariant way, 
and therefore the method would perform identically in 
non-dimensional “code units”. 

The final test, the W7 model, uses a very nontrivial 
initial condition for the ejecta that is read from a ta¬ 
ble. This test is meant to demonstrate performance on 
a pragmatic, complex problem that is not built out of 
simple analytic functions. 

3.1. Test #1: The Idealized Supernova 

This is the simplest RT-unstable supernova model that 
one can write down. A uniform-density ball of gas ex¬ 
pands freely into a uniform-density medium. Initial con¬ 
ditions are given by 

^ fl/(4rf/3) .<« 
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Fig. 4.— ID angle-averaged profiles of the idealized supernova 
at time t = 1.0. All variables are plotted for ease of comparison 
between ID, 2D and the ID model introduced in this work. The 
most dramatic multidimensional effects are the seen in the density 
and passive scalar. 


_ f ^max ij" / R') T R 

^\0 r>R 

and the pressure is negligible: 

p = io-VL.- 


(41) 


(42) 


The remaining constants are R = 0.01, Vmax = ylO/S. 
These initial conditions are chosen to give the ejecta a 
total mass M = 1 and energy E = 1, and to have it 
collide with a uniform CSM of density pcSM = 1- The 


Fig. 5.— Snapshots of the idealized supernova at various times. 
At early times {t = 0.2), the coarse-grained dynamics are only 
crudely approximated b y th e ID model. This is most likely due to 
the fact that equation II25II is used for the growth rates, whereas 
the linear growth rate would probably be more appropriate at early 
times. Nevertheless, the position of the reverse shock is accurately 
reproduced at all times, and the disruption of the extremely sharp 
contact discontinuity is reasonably well approximated. 


solution is then evaluated at time t = 1, when the ejecta 
and CSM masses are comparable, and the RT instability 
is prominent (Figure [2|). 

Figure [3] shows ID angle-averaged values of p and X as 
a function of radius. The passive scalar is plotted to show 
how much mixing is present between ejecta and CSM. In 
the ID case, there is of course no mixing. The 2D case, 
however, shows a great deal of mixing between ejecta 
and CSM, completely smoothing away the sharp density 
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Fig. 6. — Growth of the instability at early times is compared 
with analytical predictions. Given the acceleration ^ as a function 
of time, se mi- analyt ical curves are drawn for h(t) and K,{t) using 
equations l l46ll and ll48ll . respectively. These are compared with 

the numerically measured h(t) = ^.nd where 

brackets denote an average over the mixing region (equation I50II . 

jump seen in ID. The ID RT model captures this mixing 
extremely well, though again this is not too surprising, as 
the constants A, B, and C were found by trying to match 
this particular solution as closely as possible. At r = 
0.85, the blastwave is mixed with 36% ejecta according 
to the 2D results, and 37% ejecta according to the ID 
model. The model predicts slightly more penetration of 
the RT fingers (the passive scalar X is nonzero over a 
larger range in the model), but is otherwise in agreement 
with the 2D results. 

For a more complete comparison. Figure |4] plots all 
quantities as a function of radius for the ID, 2D and 
“ID+Model” calculations. The most dramatic differ¬ 
ences between ID and 2D are seen in the density p and 
the ejecta fraction V, which is why those two quanti¬ 
ties have been chosen to represent the results in most 
figures. The turbulent velocity represented by the quan¬ 
tity K is also plotted in the bottom panel of Figure IH 
This quantity is compared directly with the 2D results 
using equation (jS]) to calculate k, and evaluating the ki- 
netic turbulent energy densi ty in the same manner as in 
iDiiffell fc MacFadvenI (j2013l l. but in the nonrelativistic 
limit. Putting these results together, k is computed via 
the formula: 


the relativistic version of t his formula is presented in 
iDuffell fc MacFadveiJ (|2013ll . 

The constants A, B, and C have been tuned so that the 
profile of k matches the 2D results as closely as possible. 
These choices of A = 1.18x10“®, B = 1.2, and C = 0.102 
are fixed when applied to the remaining test problems 
(sections 13.21 - 13.51) . 

Due to the inherent scale-invariance of hydrodynamics, 
a great deal of parameter space is already implicitly ex¬ 
plored in this single test. Performing the same test with 
a different value of E, M, or pcSM would yield identi¬ 
cal results, assuming the limit R <C {M/pcsmY^^■ Since 
the ID model is introduced in a scale-invariant way, this 
model would also perform identically well for an ideal¬ 
ized supernova with any other choice of A, M, or pcSM- 
This is assuming the snapshot is taken at the re-scaled 
time 


t = to = M^/^E-^^^Pcli!. (44) 

Therefore, the only true free parameter in this test 
problem is the time at which the solution is evaluated 
(i.e. the dimensionless parameter t/to). For additional 
free parameters, one could also vary the detailed profile 
of the ejecta or CSM, but that is precisely what will be 
tested in the remaining four test problems (sections 13.21 

Figure E] presents the idealized supernova taken at sev¬ 
eral different snapshots in time, to explore how well the 
model captures the evolution of the instability, from early 
times when the density jump is very large, into late times 
when all of the gradients are smoothed out and shallow. 
At very early times {t = 0.2), the ID model only crudely 
approximates the solution, capturing the disruption of 
the extremely narrow contact discontinuity in an order- 
of-magnitude sense. Inaccuracies here may have to do 
with the fact that equation (E51) is used for the growth 
rate, which is possibly not valid at early times. 

Nevertheless, after this brief transient phase, the ID 
model rapidly approaches the correct 2D solution, cap¬ 
turing the position of the reverse shock accurately, as 
well as the disruption of the contact discontinuity. 

Finally, the initial growth phase of the instability is 
tested against analytical predictions for the evolution of 
width of the mixing layer h{t), and for K(t). Equation 
((^51) for the growth of the instability suggests that the 
width of the mixing layer obeys 

h = 2agA, (45) 

where g = P' /p is the acceleration felt by the mixing 
layer. Therefore, h{t) can be predicted to follow 

h{t) = a f dt' [ dt"2Ag{t"). (46) 

Jo Jo 


2 (T’)cons — {P)vol 

7(7 - 1 ) (-P)vol 


For K, one can assume the growth rate of the mixing 
(43) laye r is proportional to the turbulent velocity fluctuation: 


where brackets denote an average over a spherical 
shell, Qyoi denotes a volume average, while Qcons de¬ 
notes a “conservative average”, where the total mass, 
energy, and radial momentum of the shell are calculated 
and these are converted back to the primitive variables 
of density, pressure, and velocity. A full derivation of 


h oc 5u, (47) 

for some constant of proportionality. In this case, one 
can readily predict that 

K oc dt'Ag{t')^ !(?. 


(48) 
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The quantities (g), (k), (c), and the width of the mix¬ 
ing layer h are measured as averages in the mixing region, 
using the weighting function 


w(r) = A(r)(l — A(r-)), (49) 

which only has support in the mixed region. Averaged 
quantities are calculated via 


/ w(r)Q{r)dr 
f w{r)dr 


(50) 


and the mixing layer thickness is given by the standard 
deviation of w: 


h = (51) 

h{t) and the averaged are plotted as a function 

of time in Figure [SI These are plotted against computed 
values integrating equations (HSl) and (H51) . given the ac¬ 
celeration {g{t))^ as a function of time (for simplicity, 
it is assumed that A = 1 as the density jump is very 
large at early times). Growth produced by the model is 
consistent with nonlinear RT, with a = 0.2. This a is 
significantly larger than typical values, which is further 
evidence that the early-time growth is too fast (though 
the definition of h could come with an overall dimension¬ 
less coefficient). K{t) also appears consistent with (H51) . 

3.2. Test #2: The Soft Explosion 

This is another simple set-up, but with smooth initial 
conditions that generate weaker shocks. This is designed 
for comparison with other numerical implementations, 
for which shock capturing may be more challenging. The 
shocks and density jump have very different strengths 
from the first test, so that this probes a significantly 
distinct region of parameter space. Initial conditions are 
given by the following: 


P{r) = Peiir) + pcSM{r) 

(52) 

P{r) = 2.5p(r) 

(53) 

X{r) = pej(r)/p(r) 

(54) 

v{r) = (r/Ro)A(r), 

(55) 

where Rq = 0.1 and 


Pej(r)=e-(”/«°)Vi?o, 

(56) 

Pcsm(?’) = 1- 

(57) 


The initial conditions are thus completely specified by 
smooth functions of radius, and the non-negligible pres¬ 
sure will cause the shocks to be weaker and more easily 
resolved. 

Solutions for these initial conditions are presented in 
FigurejT] For this test problem, the influence of RT is not 
as significant as in the previous test, and the model cor¬ 
rectly captures this, and does not overmix. The density 
peak in ID is at about Pmax = 4.3, whereas RT instability 
reduces the density peak to Pmax ~ 3. In this case, the 
reverse shock is relatively unaffected by the turbulence. 



0 0.5 1 1.5 2 2.5 3 


Fig. 7.— ID angle-averaged profiles of the soft explosion at 
t = 1.0. Shocks are weaker and the solution is smoother than the 
previous test. RT mixing is also not as significant for this test 
problem, and this fact is faithfully recovered in the model. The 
density peak is reduced from pmax = 4.3 to pmax ~ 3. 

3.3. Test The Persistent Engine 

This test is significantly distinct from the others, as 
the energy is input to the system on a non-negligible 
timescale, modeling an engine which is persistent. This 
has the effect of producing a much more complicated 
ejecta structure, and therefore this is a much “messier” 
test than the previous ones, but more importantly the in¬ 
terior pressure causes an RT instability which is “upside- 
down” with respect to the others; there is a negative 
pressure gradient coinciding with a positive density gra¬ 
dient. In other words, low density material is accelerating 
the high-density material in front of it, in contrast with 
the low-density material decelerating high-density mate¬ 
rial behind it, as is the case in all of the other tests. In 
addition, there will be the usual RT caused by the decel¬ 
eration of the ejecta, so that the ejecta is unstable from 
both sides. This test will help to demonstrate whether 
this model is reasonable for capturing different scenarios 
which can cause RT, other than a decelerating shell of 
ejecta. 

Initial conditions are given by a static fluid of constant 
pressure: 


P = 0.1 

(58) 

u = 0 

(59) 
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Fig. 8.— ID angle-averaged profiles of the persistent engine test 
at t = 1.0. Solutions are given for two different values of the pa¬ 
rameter a, which measures the thickness of the transition between 
ejecta and CSM, setting the steepness of the density gradient in 
the unstable region. 


p= 1 + 


tanh((R — r)/a) 
2R3 


(60) 


where i? = 0.1 and a is a free parameter which controls 
the steepness of the contact discontinuity and therefore 
the growth of RT. 

Energy is deposited into the system thermally, as a 
source term. The power per unit volume is given by: 


*• = iriipys’^ ^‘ ^ 

where Re = 0.05 and T = 0.1. This results in an 
outflow with total energy E = \ when the engine has 
terminated. 

Solutions are presented in Figure [8] for two different 
values of the parameter a. For broad initial conditions 
with a = 0.05, the density gradients are shallow and RT 
is weak; the internal RT caused by the engine reduces the 
density peak by about a factor of two, and the external 
RT caused by deceleration has very little effect on the 
solution. This is reproduced in the ID model. 

The more extreme, “sharp” jump with a = 0.0005 
(lower panel of Figure [5]) results in complex, fine struc¬ 
ture in the ID solution (note: these sharp ups and downs 
are not numerical; they are resolved and fully converged 
in ID). 



8e+15 1e+16 1e+16 1e+16 2e+16 

r (cm) 


Fig. 9.— ID angle-averaged profiles of the brief encounter at 
t = 50 days. Similar to Figures 0 and 0 density is plotted along¬ 
side the passive scalar X. The multidimensional solution finds two 
distinct regions that are RT-unstable, associated with the two con¬ 
tact discontinuities in the initial condition. The RT model pre¬ 
sented here is able to capture both RT-unstable regions. 


For this value of a, RT is much stronger and diffuses 
the entire region. RT is prominent from both the inte¬ 
rior and exterior of the ejecta, and it smoothes out the 
solution so that the density peaks at around p ^ 20 in¬ 
stead of sharply peaking at p ^ 50 as in the ID version. 
It is reassuring that all of these features are faithfully 
captured in the model, as the steepness of the contact 
discontinuity is varied by several orders of magnitude. 

3.4. Test The Brief Encounter 

This is a slightly more complicated initial condition. 
The ejecta is given by a broken power-law, such that most 
of the ejected mass is contained within the radius i?i, but 
RT can be seen at early times, before the reverse shock 
passes this radius. Additionally, the encounter does not 
last forever, as this is meant to represent a finite shell of 
gas that the supernova collides with. The initial condi¬ 
tions are 


P = 


p,(Ri/i?c)-'°(r/i?i)-i 

Pc(r/i?c)“^° 

PCSM 

lO'^'pcSM 


r < Ri 
Ri < r < Rc 
Rc < r < R 2 
r > R 2 


(62) 
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Fig. 10.— Initial conditions for the W7 Model. This final test 
was read from a table, rather than being built out of simple analytic 
functions, in contrast with the previous tests. 


"10 


r < Rc 
r > Rc 


(63) 


where i?i = 4.4 x 10^"^ cm, Rc = 1.4 x 10^® cm, i ?2 = 
2.8 X 10^® cm, pcSM = 1-3 x 10“^^ g/cm^, = 4 x 10“^^ 
g/cm^, to = 5 days. For this problem, the temperature 
is negligible: 


T = 10'‘iF (64) 

The solution is evaluated at t = 45 days later, a to¬ 
tal of t -I- to = 50 days after the supernova (Figure IH]). 
At this time, the shocks have only encountered the steep 
density gradient, but the forward shock has over¬ 
taken the entirety of the shell. After crossing the shell, a 
second forward shock is pushed forward into the lower- 
density “vacuum” state. A 2D study of this solution 
shows that two separate regions are RT-unstable, associ¬ 
ated with the two different contact discontinuities in the 
initial condition. The first unstable region corresponds 
to the inner edge of the shell, being disrupted by RT 
as the shell decelerates the ejecta. The second unstable 
region applies to the interface between the outer edge 
of the shell and the “vacuum” state just outside; this 
density drop is also unstable. These two separate un¬ 
stable regions are faithfully reproduced in the ID model 
presented in this work, again using the same parameters 
calibrated by the first test. 

3.5. Test #5: The W7 Model 

W7 is a commonly-used pure deflagration type-la su¬ 
pernova ejecta model that was calcu lated using detailed 
nucleosynthesis (iNomoto et al.lHO^ . The data is pub¬ 
licly available and therefore makes an ideal test case for 
complicated supernova ejecta profiles. In this example, 
the ejecta is assumed to be expanding into a wind, with 
density 


Pwindir') — Rjv 


(65) 



1.2e+20 1.5e+20 1.8e+20 2.1e+20 

r (cm) 


Fig. 11.— W7 Model at time t = 1187 years. At this late time, 
the supernova remnant has swept up a signficant portion of the 
CSM. 

Results from the W7 model are shown in Figure [TTJ 
The supernova remnant is shown at very late times, 
t = 1187 years, after which a substantial amount of the 
flow has been decelerated by the wind (it might be more 
realistic to use a shallower density profile at such large 
radii, but this is an unimportant detail for the purposes 
of the current study). The ID solution is characterized 
by a sharp density jump that is unstable in 2D. Again, 
the model presented in this work accurately captures the 
2D effects, and provides a much more accurate approxi¬ 
mation to the 2D solution than the ID version. 

In the ID case, there is again a sharp contact discon¬ 
tinuity, but in 2D, there is significant mixing between 
ejecta and CSM. This also pushes the reverse shock sig¬ 
nificantly backward. All of these basic features are repro¬ 
duced by the ID model, again using the same parameters 
A, B, and C which were calibrated by the first test. The 
main distinction is again the penetration depth of the RT 
fingers. At radii up to r ~ 1.7 x 10^° cm, the model pre¬ 
dicts the correct amount of mixing (about 10% of ejecta 
mixed with CSM), but at larger radii the model some¬ 
what overpredicts how deeply the fingers penetrate when 
compared to the 2D results. 


with AT = 5 X 10^^ g/cm. This wind corresponds to a 
mass loss rate of per year and a wind velocity 

of 1000 km/s. The initial density and velocity profile is 
plotted in Figure ITOl 


4. SUMMARY 

Connecting supernova models to observations may 
require extensive searches of parameter space, which 
may only be possible in ID, especially if these su- 
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Fig. 12.— Test #1 is performed at five different resolutions in 
2D, in order to test convergence. In all cases, 59 ~ 5r/r. 

pernova models include full radiation hydrodynamics, 
which is expensive in multidimensions. The model pre¬ 
sented here provides a useful tool for exploring this 
parameter space. This is potentially a big improve¬ 
ment over the current standard, which is to perform 
the calculation in ID, then man ually smooth out offend¬ 
ing regions in post-processing (IPinto fc WooslevI 119881 : 
iKasen fc Woosle\^l2009t iHeger fc Wooslevll2010ll . 

Since the modihed equations give a value for k = 
(5u^/Cs, this potentially provides an estimate for the 
magnetic field strength, if magnetic fields are as¬ 
sumed to be in equipartition with kinetic fluctuations 


(iHaueen et al.ll2003l ISchekochihin et al.ll2004 IBeresnvak 

2012 

; Zrake & MacFadvenI 120131 iDuffell & MacFadven 

2013 

.120141). In this case, er ~ (t(t — 1)/2)k, by Equa- 


tion ([8]). 


The current model is accurate enough to be used im¬ 
mediately in ID studies of supernovae, but it should be 
stressed th at th is is on ly a first step (actually a second 
step, after iGullI (IlDTSlD toward a complete RT model. 
It should also be noted that the model includes several 
fitting parameters, and that this might not be a unique 
method of matching the multidimensional results. 

The model is designed and calibrated to match late¬ 
time behavior, after initial transients have washed out 
and the swept-up mass is comparable to the mass of the 
ejecta. However, at very late times, the decay of turbu¬ 
lence is important, and this has not yet been accounted 
for. A more accurate model would utilize a linear model 
for the early-time growth rate, and include the decay 
term to match late-time evolution. A connection between 
the early and late phas es of growth has been realized in 
the ID ODE models of (lRanishawlll998ll . and this model 
was incorporated in a multi-phase PDE mixing model 
(|R,amshawl l200fl[ l . so it is likely possible to incorporate 
such modifications into the single-fluid PDE model pre¬ 
sented here in a future study. 

Another important point is that the constants A, B, C 
were calibrated using 2D numerical results. 3D results 
would be much more accurate, even if they are not re¬ 
solved as well as their 2D counterparts. It would also be 
desirable to extend this model to relativistic flows, for 
application to gamma ray burst afterglows. Finally, it is 
well-known that cooling from cosmic rays can strongly 
impact the dynamics of RT, and it is unclear whether 
this would need to be accounted for in the model. This 
should be tested explicitly in a future study. 
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the NASA Advanced Supercomputing (NAS) Division 
at Ames Research Center. I am grateful to Chris Mc¬ 
Kee, Dan Kasen, Eliot Quataert, Andrei Cruzinov, and 
Nathan Roth for helpful comments and discussions. I 
would especially like to thank Chelsea Harris for encour¬ 
aging me to attempt this study in the hrst place, and Bill 
Paxton for additional encouragement and for inspiring 
several of the tests performed in this study. I would also 
like to thank the anonymous referee for the extremely 
detailed and thorough review. 

APPENDIX 

RESOLUTION STUDY 

In order to demonstrate that the tests in this study 
are based on converged numerical solutions, a resolution 
study is performed o n th e first test (the “Idealized Su¬ 
pernova” , see section IQ and Figure [2|). 

This test is performed in 2D at five different resolu¬ 
tions, shown in Figure 1121 Even the lowest resolution 
Ar/r ~ 0.02 captures the coarse-grained effects reason¬ 
ably well (surprisingly so). As the resolution is increased, 
the angle-averaged quantities measured in this study ap¬ 
proach a converged solution. 
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